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ABSTRACT 



>, 
C^ I We present a systematic treatment of the linear theory of scalar gravitational 

perturbations in the synchronous gauge and the conformal Newtonian (or longitudinal) 

gauge. We first derive the transformation law relating the two gauges. We then 

write down in parallel in both gauges the coupled, linearized Boltzmann, Einstein and 

0^ ■ fluid equations that govern the evolution of the metric perturbations and the density 

fluctuations of the particle species. The particle species considered include cold dark 

matter (CDM), baryons, photons, massless neutrinos, and massive neutrinos (a hot 

2 ! dark matter or HDM candidate), where the CDM and baryon components are treated 

^ I as fluids while a detailed phase-space description is given to the photons and neutrinos. 

The linear evolution equations presented are applicable to any = 1 model with CDM 

or a mixture of CDM and HDM. Isentropic initial conditions on super-horizon scales 

^ ' are derived. The equations are solved numerically in both gauges for a CDM+HDM 

j^ ■ model with f^coid = 0.65, ilhot = 0.3, and ribaryon = 0.05. We discuss the evolution 

of the metric and the density perturbations and compare their different behaviors 

outside the horizon in the two gauges. In a companion paper we integrate the geodesic 

equations for the neutrino particles in the perturbed conformal Newtonian background 

metric computed here. The purpose is to obtain an accurate sampling of the neutrino 

phase space for the HDM initial conditions in A^-body simulations of the CDM+HDM 

models. 
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1. Introduction 



The theory of galaxy formation based on gravitational instability is aimed at describing how 
primordially-generated fluctuations in matter and radiation grow into galaxies and clusters of 
galaxies due to self-gravity. A perturbation theory can be formulated when the amplitudes of the 
fluctuations are small, and the growth of the fluctuations can be solved from the linear theory. 
Such linear theory for a perturbed Friedmann-Robertson- Walker universe was first developed by 
Lifshitz (1946), later reviewed in Lifshitz & Khalatnikov (1963). The subsequent work can be 
found summarized in the textbooks by Weinberg (1972) and Peebles (1980), and in the reviews by 
Kodama & Sasaki (1984) and Mukhanov, Feldman & Brandenberger (1992). 

In the early universe, gravitational perturbations are inflated to wavelengths beyond the 
horizon at the end of the inflationary epoch. Fluctuations of a given length scale reenter the 
horizon at a later time when the horizon has grown to the size of the fluctuations. Although the 
process of galaxy formation in recent epochs is well described by Newtonian gravity (and other 
microphysical processes such as hydrodynamics), a general relativistic treatment is required for 
perturbations on scales larger than the horizon size before the horizon crossing time. The use of 
general relativity brought in the issue of gauge freedom which has caused some confusion over 
the years. Lifshitz (1946) adopted the "synchronous gauge" for his coordinate system, which has 
since become the most commonly used gauge for cosmological perturbation theories. However, 
some complications associated with this gauge such as the appearance of coordinate singularities 
and spurious gauge modes prompted Bardeen (1980) and others (e.g. Kodama & Sasaki, 1984) 
to formulate alternative approaches that deal only with gauge-invariant quantities. A thorough 
review of the gauge-invariant perturbation theory and its appliaiton to texture-seeded structure 
formation models is given by Durrer (1993). Another possibility is to adopt a different gauge. We 
will discus in detail in this paper the conformal Newtonian (or the longitudinal) gauge (Mukhanov 
et al. 1992), which is a particularly convenient gauge to use for scalar perturbations. 

This paper serves two purposes. First, it is an independent paper in which we present and 
compare a systematic treatment of the linear theory of scalar isentropic gravitational perturbations 
in the synchronous and conformal Newtonian gauges. The coupled, linearized Einstein, Boltzmann, 
and fluid equations for the metric and density perturbations are presented in parallel in the two 
gauges. We give a full discussion of the evolution of five particle species: cold dark matter (CDM), 
hot dark matter (HDM, i.e., massive neutrinos), baryons, photons, and massless neutrinos. The 
CDM and the baryon components behave like collisionless and collisional fluids, respectively, 
while the photons and the neutrinos require a phase-space description governed by the Boltzmann 
transport equation. We also derive analytically the time dependence of the perturbations on 
scales larger than the horizon to illustrate the dependence on the gauge choice. This information 
is needed in the initial conditions for the numerical integration of the evolution equations. The 
linear theory discussed in this paper can be applied to il = 1 models with various dark matter 
compositions, e.g., pure CDM, pure HDM, or a mixture of CDM and HDM. 



This paper also serves as a companion paper to Ma & Bertschinger (1994) in which we 
reported the main results from our linear calculation of the full neutrino phase space in a 
CDM+HDM model with Qcoid = 0.65, Ohot = 0.3, l^baryon = 0.05, and Hq = 50 km s"^ Mpc^^. 
(The corresponding neutrino mass is rrii, = 6.985 eV.) The motivation was to obtain an accurate 
sampling of the neutrino phase space for the HDM initial conditions in A^-body simulations of 
structure formation in CDM+HDM models. We adopted a two-step Monte Carlo procedure to 
achieve this goal: (1) Integrate the coupled, linearized Boltzmann, Einstein, and fluid equations 
for all particle species in the model (i.e., CDM, HDM, photons, baryons and massless neutrinos) 
to obtain the evolution of the metric perturbations; (2) Follow the trajectories of individual 
neutrinos by integrating the geodesic equations using the metric computed in (1). Since no 
coordinate singularities occur in the conformal Newtonian gauge and the geodesic equations have 
simple forms, the geodesic integration in step (2) was carried out in this gauge, starting shortly 
after neutrino decoupling at redshift z ~ 10^ until z = 13.5. We focus on step (1) in this paper. 
Following historical precedents, we first developed the code for the Boltzmann integration in 
the synchronous gauge. The transformation relating the synchronous gauge and the conformal 
Newtonian gauge was then derived and used to compute the metric perturbations in the latter 
gauge for step (2) of the calculation. Subsequently we developed a code to perform the full 
integration in the conformal Newtonian gauge. 

The organization of this paper is as follows. In §2 we write down the metric for the two gauges 
and summarize their properties. In §3, we derive the gauge transformation relating two arbitrary 
gauges and obtain the transformation between the synchronous and the conformal Newtonian 
gauges. The linearized evolution equations for the metric and the density perturbations are given 
in §§4 and 5. Section 4 discusses the Einstein equations with emphasis on the source terms, the 
energy-momentum tensor, in the two gauges. The perturbed fluid equations are derived from the 
energy-momentum conservation, which are applied to CDM and the baryons in §5. The rest of §5 
contains detailed treatments of the photon and neutrino phase space distributions, recombination, 
and the coupling of photons and baryons. The photon and neutrino distribution functions are 
expanded in Legendre polynomials, reducing the linearized Boltzmann equation to a set of coupled 
ordinary differential equations for the expansion modes. The massive neutrinos require a slightly 
more complicated treatment due to the nontrivial time dependence of the energy-momentum 
relation. Section 6 discusses the behavior of the perturbations before horizon crossing. The 
necessary initial conditions for the variables in the two gauges are given. Section 7 presents the 
numerical results for the evolution of the perturbations in our CDM+HDM model. 



2. The Tv^^o Gauges 

We consider only spatially flat (Q = 1) background spacetimes with isentropic scalar metric 
perturbations. The spacetime coordinates are denoted by x'^, jx G (0,1,2,3), where x^ is the 



time component and x*, i £ (1,2,3) are the spatial components in Cartesian coordinates. Greek 
letters a, (3, 7 and so on always run from to 3, labeling the four spacetime-coordinates; Roman 
letters such as i,j, k always run from 1 to 3, labeling the spatial parts of a four-vector. Repeated 
indices are summed. Since our interests lie in the physics in an expanding universe, we use 
comoving coordinates x^ = (r, x) with the expansion factor a{T) of the universe factored out. The 
comoving coordinates are related to the proper time and positions t and r by dx^ = dr = dt/a(T), 
dx = dr/a{T). Dots will denote derivatives with respect to r: d = da /dr. The speed of light c is 
set to unity. 

The components qqq and gQi of the metric tensor in the synchronous gauge is by definition 
unperturbed. The line element is given by 

ds^ = a^{T){-dT'^ + {5ij + hij)dx'dx^} . (1) 

The metric perturbation hij can be decomposed into a trace part h = ha and a traceless part 
consisting of three pieces, /i| •, hjj, and hf,, where hij = h5ij/3 + hj- + hj-, + hj,. By definition, the 
divergences of h'j- and h^ (which are vectors) are longitudinal and transverse, respectively, and hj- 
is transverse, satisfying 

eijkdjdihjf^ = , didjhjj = , dihjj = . (2) 

1 be written in term 
divergenceless vector A as 



It then follows that /i!' • can be written in terms of some scalar field A and hj-^ in terms of some 



4 = {d^d.-U^.V'^A, 

hij = diAj + djAi, d,Ai = 0. (3) 

The two scalar fields h and A (or h--) characterize the scalar mode of the metric perturbations, 
while Ai (or hjj) and hfj represent the vector and the tensor modes, respectively. 

We will be working in the Fourier space k in this paper. We introduce two fields h{k, r) and 
r]{k,T) in /c-space and write the scalar mode of hij as a Fourier integral 

hij{x, t) = d^ke"-^"-^ \ kikjh{k, r) + {kikj 6ij) 6ri{k, t)\ , k = kk. (4) 

Note that h is used to denote the trace of hij in both the real space and the Fourier space. 

In spite of its wide-spread use, there are serious disadvantages associated with the synchronous 
gauge. Since the choice of the initial hypersurface and its coordinate assignments are arbitrary, 
the synchronous gauge conditions do not fix the gauge degrees of freedom completely. Such 
residual gauge freedom is manifested in the spurious gauge modes contained in the solutions to 
the equations for the density perturbations. The appearance of these modes has caused some 
confusion over the years and prompted Bardeen (1980) to formulate alternative approaches that 



deal only with gauge- invariant quantities. Another difficulty with the synchronous gauge is that 
since the coordinates are defined by freely falling observers, coordinate singularities arise when two 
observers' trajectories intersect each other: a point in spacetime will have two coordinate labels. 
A different initial hypersurface of constant time has to be chosen to remove these singularities. 

The conformal Newtonian gauge (also known as the longitudinal gauge) advocated by 
Mukhanov et al. (1992) is a particularly simple gauge to use for the scalar mode of metric 
perturbations. The perturbations are characterized by two scalar potentials ■0 and (p which appear 
in the line element as 

ds^ = a'^{T) {-(1 + 2-0)^x2 + (1 - 2<P)dx'dxi^ . (5) 

It should be emphasized that the conformal Newtonian gauge is a restricted gauge since the metric 
is applicable only for the scalar mode of the metric perturbations; the vector and the tensor 
degrees of freedom are eliminated from the beginning. Nonetheless, it can be easily generalized 
to include the vector and the tensor modes. We will confine our discussion here to the scalar 
perturbations only. 

One advantage of working in this gauge is that the metric tensor g^^j is diagonal. This 
simplifies the calculations and leads to simple geodesic equations (Ma & Bertschinger 1994). 
Another advantage is that ip plays the role of the gravitational potential in the Newtonian limit 
and thus has a simple physical interpretation. Moreover, the two scalar potentials ip and (j) in 
this gauge are identical (up to a minus sign) to the gauge-invariant variables $a and ^h Bardeen 
constructed (1980). No gauge modes are present in this gauge to obscure the meaning of the 
physical modes since the gauge freedom is completely fixed for Vt = 1 aside from the addition of 
spatial constants to ip and (j). The second scalar (j) is required when the energy-momentum tensor 
T^^i, contains a nonvanishing traceless and longitudinal component. As we will see in equation 



(21) and the Einstein equation ( |22d[ ), this component provides the source term for the constraint 
equation for {(j) — tp). When this component is absent, the two scalar potentials ijj and (p are 
identical. 



3. Gauge Transformations 

In this section we first derive the transformation law relating two arbitrary gauges. From it, 
the gauge transformation relating the synchronous gauge and the conformal Newtonian gauge is 
readily obtained. 

A perturbed flat Priedmann-Robertson- Walker metric can be written in general as 
500 = -a\T){l + 2^P{x,T)} , 

90i = -a^{T)wi{x,T), (6) 

Qij = a^{T) {[1 - 2(l){x,T)]5ij + (Tij(f,r)} , CTm = 



where the functions ip, (j), Wi and cjij represent metric perturbations about the Robertson- Walker 
spacetime and are assumed to be small compared with unity. The trace part of the perturbation 
to Qij is absorbed in 0, and aij is taken to be traceless. 

Consider a general coordinate transformation from a coordinate system x'^ to another x^ 

x'' ^x^' = x^' + d^'{x''). (7) 

We write the time and the spatial parts separately as 

X = X + a{x,T) , 
^ = x + V(3{x,T)+e{x,T), V-€ = 0, (8) 

where the vector d has been decomposed into a longitudinal component V/3 (V x V/3 = 0) and a 
transverse component e* (V • e*= 0). The requirement that ds^ be invariant under this coordinate 
transformation leads to 

9fMu{x) = g^iv{x) - g^f3{x)dyd'^ - gau{x)d^d" - d"dagfMu{x) + Oid^) . (9) 

We note that both sides of this equation are evaluated at the same coordinate values x in the two 
gauges, which do not correspond to the same physical point in general. Assuming d^ to be of 
the same order as the metric perturbations 'ip,Wi,(p ^-nd Uj^, the metric perturbations in the two 
coordinate systems are related to first order in the perturbed quantities by 

^(x, r) = ip{x,T) — a{x,T) a{x,T), (10a) 

a 

Wi{x,T) = Wi{x,T) +dia{x,T) - di${x,T) - ei{x,T) , (10b) 

0(x,r) = 0(x,r) + -V2/3(x,r) + -a(x,r), (10c) 

6 a 

a,j{x, r) = aij{x, r) - 2 | (diOj - ^^ij^^) P{x, r) + ^ (d^e^ + O^eA . (lOd) 

We can further decompose the transformations of Wi and aij above into longitudinal and 
transverse parts: 

wj {x, t) = w] {x, t) + dia{x, r) - diP{x, r) , 

wH^^'t) = w^{x,T)-ei{x,T), (11) 



and 
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(f,r) = al{x,T)-2(^d^d,-^6^jV^^p{x,T). 



^ijix,r) = alj{x,T)-{di €j + dj Ci ) , 

^ij{x,r) = afj{x,T), (12) 



where Wi = w^ + w^- , atj = o"] • + a^j + afj , and (t| ■ , a^ and afj obey equation (^) . Equations 
([ToD-(p!2[) describe the transformation of metric perturbations under a general infinitesimal 
coordinate transformation. 



We can now use equations (|10| ) to relate the scalar metric perturbations {(j), ip) in the 
conformal Newtonian gauge to hij = h6ij/3 + hj- in the synchronous gauge. Let x'^ denote the 
synchronous coordinates and x^ the conformal Newtonian coordinates with x^ = x^ + d^. Prom 
equations (pA|) and (|T^), we find 



a{x,T) = /3(x, r) + x(t) , 

ei{x,T) = ei{x) , 

hlj{x,T) = '^(^^ 1;:.V72 

diej + dj€i = 0, 



2{didj--6ijV']P{x,T) 



(13a) 
(13b) 

(13c) 

(13d) 



where xi^) is an arbitrary function of time, refiecting the gauge freedom associated with the 
coordinate transformation: x^ = x^ + x(''")) ^* = x*. This transformation corresponds to a global 
redefinition of the units of time with no physical significance; therefore we shall set x = from 
now on. From equations ( |10a| ) and ( |10c| ) we then obtain 



(14) 



^p{x,T) 
0(x,r) 



+/3(x,r) + -/3(x,r), 

-i/i(x,r)- V/?(x,r)--/?(f,r) 
DO a 



and P is determined by h" in equation ( |13c) ). 

In terms of h and 77 introduced in equation (^), /i^ • in the synchronous gauge is given by 

hl{x, t)= j (fk e^^-^ (kk - -6ij) [h{k, r) + 67]{k, r)} , k = kk . 



(15) 



Comparing /i"- in equations (|13c|) and (|15|), we can read off /?: 

/3(x,r) =y'd3^e*^^-^{/i(^,r) + 6r/(fe,r)} . (16) 

Then from equations (14), the conformal Newtonian potentials (p and ip are related to the 



synchronous potentials h and r/ in /c-space by 



V'(^,t) = ^|/i(^,r)+6f?(^,r) + ^[/i(fc,r)+677(^,r) 

,r) = r]{k,T)-—-r- h{k,T)+6r]{k,T) 
Ik"^ a L 



(17) 



The other components of the metric perturbations, Wi, ajj , and ajj , are zero in both gauges. 



4. Einstein Equations and Energy-Momentum Conservation 

For a homogeneous Friedmann-Robertson- Walker universe with energy density p{t) and 
pressure P{t), the Einstein equations give the following evolution equations for the expansion 
factor a{T): 

-Y = ^GaV--K, (18) 

a/ 6 



-Ga\-p + W), (19) 

dr \aj 6 

where the dots denote derivatives with respect to r, and k is positive, zero, or negative for a 
closed, flat, or open universe, respectively. We consider only = 1 models in this paper, so we set 
K = 0. It follows that the expansion factor scales as a oc r in the radiation-dominated era and 
a oc r^ in the matter-dominated era. 

We find it most convenient to solve the linearized Einstein equations in the two gauges in 
the Fourier space k. In the synchronous gauge, the scalar perturbations are characterized by 
h{k,T) and ri{k,T) in equation (Q). In terms of h and r], the time-time, longitudinal time-space, 
trace space-space, and longitudinal traceless space-space parts of the Einstein equations give the 
following four equations to linear order in /c-space: 



Synchronous gauge 



k'^r]---h = ATrGa^6T%{Syn), (20a) 
2 a 

k^ = 4TTGa^{p + P)e{Syn), (20b) 

h + 2-h-2k^r] = -8TTGa^6T\(Syn), (20c) 
a 

h + 6ii + 2-(h + 6^) -2k^r] = -2A-KGa^{p + P)e{Syn) . (20d) 



a 

The label "Syn" is used to distinguish the components of the energy-momentum tensor in the 
synchronous gauge from those in the conformal Newtonian gauge. The variables 9 and G are 
defined as 

(p + p)0 = ik^T^^j , {-p + P)e = -{kk, - -5ij)T.'j , (21) 

and S j = T j — 5'^jT^k/'^ denotes the traceless component of Tj. When the different components 
of matter and radiation (i.e., CDM, HDM, baryons, photons, and massless neutrinos) are treated 
separately, {p + P)6 = J2iiPi + Pi)(^i and {p + P)@ = J2i{Pi + Pi)&i j where the index i runs over 
the particle species. 

In the conformal Newtonian gauge, the first-order perturbed Einstein equations give 
Conformal Newtonian gauge — 

k'^(l) + 3-(4> + -'ili) = 4TTGa^6T%{Con) , (22a) 

a \ a J 





\ a J 


= A7rGa\p + P)e{Con) , 


(22b) 


-('0 + 2(/))+ (2-- 
a \ a 


-— V' + ^(</'-V') = 

(X j o 


= yGa^^T^Con), 


(22c) 




e{^-^p) -- 


= 12TTGa^{p + P)e{Con), 


(22d) 



where "Con" labels the conformal Newtonian coordinates. Next we will derive the transformation 
relating 6T^i, in the two gauges. 

For a perfect fluid of energy density p and pressure P, the energy-momentum tensor has the 
form 

T^ = Pg''^ + {p + P)U''Uu, (23) 



where U^ = dx^/V—ds"^ is the four-velocity of the fluid. The pressure P and energy density p of 
a perfect fluid at a given point are defined to be the pressure and energy density measured by a 
comoving observer at rest with the fluid at the instant of measurements. For a fluid moving with 
a small coordinate velocity v'' = dx'^/dr, v^ can be treated as a perturbation of the same order as 
6p = p — p, 5P = P — P, and the metric perturbations. Then to linear order in the perturbations 
the energy-momentum tensor is given by 



T% -- 


= -{p + 6p), 


T% - 


= {p + p)v, = -r 


Tj = 


= iP + 5P)5'j + ^ 



1 

i 



0, (24) 



where we have allowed an anisotropic shear perturbation S j in T*,-. As we shall see, since the 
photons are tightly coupled to the baryons before recombination, the dominant contribution to 
this shear stress comes from the neutrinos. We note that for a fluid, 9 defined in equation ( |2l| ) is 
simply the divergence of the fluid velocity: 9 = ik^Vj. 

The energy-momentum tensor T^iy{Syn) in the synchronous gauge is related to T'';^(Con) in 
the conformal Newtonian gauge by the transformation 

where x^ and x^ denote the synchronous and the conformal Newtonian coordinates, respectively. 
It follows that to linear order, rOo(Syn) = r°o(Con) ,r°j(Syn) = r°j(Con) + ikja{p + P) , and 



T j(Syn) = T j(Con) , where a = x^ — x'^ = {h + 6fi)/2k'^ in fc-space from equations ( [I3aD and (16). 
Let 5 = Sp/p = —6T%/p. Evaluating the perturbations at the same spacetime coordinate values, 
we obtain 

(26a) 

(26b) 
(26c) 
(26d) 



6{Sjn) 


p 


9{Syn) 


= 9{Con)-ak'^, 


SPiSyn) 


= 6P{Con)-aP, 


e(Syn) 


= e(Con). 
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This transformation also applies to individual species when more than one particle species 
contributes to the energy-momentum tensor, provided that the appropriate p and P are used for 
each component. 

The non-relativistic fluid description is appropriate for the CDM and the baryon components. 
The photon and the neutrino components, however, can be appropriately described only by their 
full distribution functions in phase space. The energy-momentum tensor in this case is expressed 
through integrals over momenta of the distribution functions. We will discuss it in detail in §5. 

The conservation of energy-momentum is a consequence of the Einstein equations. Let 
w = P/p describe the equation of state. Then the perturbed part of energy-momentum 
conservation equations 

T^"";,, = d^Tf"" + V^pT""^ + r^T'^^ = (27) 

in fc-space implies 

Synchronous gauge — 

e = -«(i-3«,)e-^i^e + ^^*t=j-t=e, (28) 

a 1 + If 1 + w 



Conformal Newtonian gauge — 

9 = -^{i-?,w)e-^^e+^-^I^kH-k'Q + k^i,. (29) 

a 1+w 1+w 

These equations are valid for a single uncoupled fluid, or for the net (mass-averaged) 5 and 9 
for all fluids. They need to be modified for individual components if the components interact 
with each other. An example is the baryonic fluid in our model, which couples to the photons 
before recombination via Thomson scattering. In the next section we will show that an extra term 
representing momentum transfer between the two components needs to be added to the 5 equation 
for the baryons. 

For the isentropic primordial perturbations considered in this paper, the equations above 
simplify since 5P/5p = c^ = w, where Cg is the adiabatic sound speed in the fluid. Although 
entropy is generated from the coupling of baryons and photons before recombination, it is a 
first-order perturbation and only appears in the energy-momentum conservation equations in the 
second order. Thus the terms proportional to {SP/6p — w) above can be neglected when the 
equations are applied to CDM and baryons in the next section. Even in the case of isocurvature 
models, which may have large entropy perturbations ab initio, these terms are generally small. 



- 11 - 

5. Evolution Equations for Matter and Radiation 

5.1. Phase Space and the Boltzmann Equation 

A phase space is described by six variables: three positions x* and their conjugate momenta 
Pi. Our treatment of phase space is based on the time-slicing of a definite gauge (synchronous or 
conformal Newtonian). Although this approach is not manifestly covariant, it yields correct results 
provided we convert gauge-dependent quantities to observables at the end of our computation. 

The conjugate momentum has the property that it is simply the spatial part of the 



4-momentum with lower indices, i.e., for a particle of mass m, Pi = mUi , where Ui = dxi/V—ds'^. 
One can verify that the conjugate momentum is related to the proper momentum p* = pi measured 
by an observer at a fixed spatial coordinate value by 

Pi = a{5ij + ^hij)p> , in synchronous gauge , 

Pi = a(l — 4>)pi , in conformal Newtonian gauge . (30) 

In the absence of metric perturbations, Hamilton's equations imply that the conjugate momenta 
are constant, so the proper momenta redshift as a^^. 

The phase space distribution of the particles gives the number of particles in a differential 
volume dx^dx^dx'^dPidP2dP2, in phase space: 

f{x\ Pj,T)dx^dx^dx^dPidP2dP3 = dN . (31) 

Importantly, / is a scalar and is invariant under canonical transformations. The zeroth-order phase 
space distribution is the Fermi-Dirac distribution for fermions (-|- sign) and the Bose-Einstein 
distribution for bosons (— sign): 

where e = a(p^ -|- mPY'"^ = (P^ -|- a'^m'^)^''^ , Tq = aT denotes the temperature of the particles 
today, the factor Qs is the number of spin degrees of freedom, and hp and k^ are the Planck and 
the Boltzmann constants. 

When the spacetime is perturbed, x* and Pi remain canonically conjugate variables, with 
equations of motion given by Hamilton's equations (Bertschinger 1993). However, following 
common practice (e.g.. Bond & Szalay 1983) we shall find it convenient to replace Pj by Qj = apj 
in order to eliminate the metric perturbations from the definition of the momenta. Moreover, 
we shall write the comoving 3-momentum Qj in terms of its magnitude and direction: Qj = qrij 
where n'^rii = dijU^n^ = 1. Thus, we change our phase space variables, replacing f{x^,Pj,T) by 
/(x*, q, nj,T). While this is not a canonical transformation (i.e., Qi is not the momentum conjugate 
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to x^), it is perfectly valid provided that we correctly transform the momenta in Hamilton's 
equations. Note that we do not transform /. Because qj are not the conjugate momenta, d^xd^q 
is not the phase space volume element, and fd^xd^q is not the particle number. In the conformal 
Newtonian gauge, for example, (1 — ?>(j)) f d^ xd^ q is the particle number; this result is sensible 
because a(l — 4))dx'^ is the proper distance. 

In the perturbed case we shall continue to define e as o(r) times the proper energy measured 
by a comoving observer, e = {q^ + a^m?Y''^ . This is related to the time component of the 
4-momentum by Pq = — e in the synchronous gauge and Pq = —(1 + ip)e in the conformal 
Newtonian gauge. For the CDM+HDM model we are interested in, the photons, the massless 
neutrinos, and the 7 eV neutrinos at the time of neutrino decoupling are all ultra-relativistic 
particles, so e in the unperturbed Fermi-Dirac and Bose-Einstein distributions can be simply 
replaced by the new variable q. 

The general expression for the energy-momentum tensor written in terms of the distribution 
function and the 4-momentum components is given by 

T^. = JdPidP2dP3 i-g)-^/^ ^f{x\Pj,T) , (33) 

where g denotes the determinant of g^u . It is convenient to write the phase space distribution as a 
zeroth-order distribution plus a perturbed piece in the new variables q and Hj: 

f{x\Pj,T)=fo{q)[l + ^{x\q,nj,T)] . (34) 

In the synchronous gauge, we have (—5)"^'^ = a~^(l — ^h) and dPidP2dPs = 
(1 + ^hijninj)q^dqdQ. to linear order, where h = ha and dil is the solid angle associated 
with direction rij. Using the relations J dQuiUj = 47r5jj/3 and J dQrii = J dQuinjUk = 0, it then 
follows from equation (|3^) that 

T% = -a-* I q^dq dQ ^ q^ + m^a^ Uq) (1 + ^) , 

T\ = a-^ Jq'^dqdnqmfoiq)^, (35) 

T) = a-' f q^dqdn^^^i=fo{q){l + ^) 

to linear order in the perturbations. Note that we have eliminated the explicit dependence on the 
metric perturbations in equation (|3^ ) by redefining Pi in terms of q and nj. Note also that the 
comoving energy e{q,T) = [q^ + a^m^)^" is used in the integrands but not in the argument of the 
unperturbed distribution function. 

In the conformal Newtonian gauge, (— 5)^"*^'^ = a~^(l — ip + 3(/>) and dPidP2dP^ = 
(1 — 3(p)q'^dqdQ. It then follows that the components of the energy-momentum tensor have the 
same form as in equations (^). Of course, it is understood that the variables q and rii in this case 
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are defined in relation to the conjugate momentum Pi in the conformal Newtonian coordinates, and 
not the synchronous coordinates. (They diff'er because comoving observers in the two coordinate 
systems are not the same.) The expansion factor a and ^ are evaluated at the coordinates (x*,r) 
in the conformal Newtonian gauge. 

The phase space distribution evolves according to the Boltzmann equation. In terms of our 
variables {x'^,q,nj,T) this is 

dr dr dr 9x' dr dq dr drii \dT ) q 

where the right-hand side involves terms due to collisions, whose form depends on the type of 
particle interactions involved. From the geodesic equation 

P°^ + r'^„«P"P'3 = 0, (37) 

dr 

it is straightforward to show that 

dq/dr = —-qhijUiUj in synchronous gauge , 

dq/dr = qcf) — e(g, r) Uidiip in conformal Newtonian gauge , (38) 

and dui/dr is 0{h). Since df/drii is also a first-order quantity, the term {dni/dT){df /drii) in the 
Boltzmann equation can be neglected to first order. Then the Boltzmann equation in fc-space can 
be written as 



Synchronous gauge 

\I/ n ^ Win f^ r /) 4- fir, „ 

-{k ■ n 



— \-i-{k- n)w + 



dr e dlnq 



h + Qfj^f ^^2 



Conformal Newtonian gauge 



-— + i-{k-n)^ + -— 

OT e dmq 



(f) — i —{k ■ n) tp 



1 (^^ 



/o \dTJc 



(39) 



(40) 



Equations ( |39[ ) and (|^) can also be derived using the canonical phase space variables re* and Pj 
and Hamilton's equations (instead of the geodesic equation) , followed by a transformation from Pj 
to qrij. 

The terms in the Boltzmann equation depend on the direction of the momentum n only 
through its angle with k. (We shall see that this is true of the collision term for photons as well 
as the convective and metric perturbation terms.) Therefore, if the momentum-dependence of the 
initial phase space perturbation is axially symmetric about /c, it will remain axially symmetric. 
If axially- asymmetric perturbations in the neutrinos or other collisionless particles are produced, 
they would generate no scalar metric perturbations and thus would have no effect on other species. 
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Therefore, we shall assume that the initial momentum-dependence is axially symmetric so that 
^ depends on q = qh only through q and k ■ h. This assumption, which effectively reduces the 
dimensionality of phase space perturbations by one (after Fourier transforming on the spatial 
coordinates), has been made (implicitly, if not explicitly) in all previous studies of the evolution of 
scalar perturbations. 



5.2. Cold Dark Matter 

CDM interacts with other particles only through gravity and can be treated as a pressureless 
perfect fluid. The CDM particles can be used to define the synchronous coordinates and therefore 
have zero peculiar velocities in this gauge. Setting 9 = Q = and w = w = c^ = in equations 
(M) leads to 



Synchronous gauge — 

5, = -h. (41) 

The CDM fluid velocity in the conformal Newtonian gauge, however, is not zero in general. In 



fe-space, equations (29) give 



Conformal Newtonian gauge 



6, = -9, + 3(^, ec = —9c + k^iJ. (42) 

a 



The subscript c in 5c and 9c denotes the cold dark matter. 



5.3. Massless Neutrinos 

The energy density and the pressure for massless neutrinos (labeled by subscripts v) are 
p,y = 3P,y = —T^o = T\. From equations (35) the unperturbed energy density p,^ and pressure P^ 



are given by 

p, = 3P, = a-4 J q^dqdn qfo{q) , (43) 

and the perturbations of energy density 5pu, pressure SP^, energy flux (JT^j , and shear stress 
T,'lj = Tlj - PySij are given by 

5p^ = 35P^ = a-^ Jq^dqdnqfoiq)^, 
8T^,i = a-^ fq^dqdnqmUq)^, (44) 



T,lj = a "^ q'^dqdQqiniUj - -6ij)fo{q)'^ 
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The unperturbed energy flux and shear stress are zero. 

The Boltzmann equation siniphfies for massless particles, for which e = q. To reduce the 
number of variables we integrate out the g-dependence in the neutrino distribution function and 
expand the angular dependence of the perturbation in a series of Legendre polynomials Pi{k ■ n): 



jq^dqqfoiq) ^ 



n] 



(45) 



«=o 



As noted in §5.1, the dependence on n arises only through k ■ h, so that a general distribution may 
be represented as in equation (^5|). 

In terms of the new variable Fy{k,n,T) and its harmonic expansion coefficients, the 



perturbations 5u, Pu, 6u, and Qu (defined in eq. ^T]) take the form 

5y = — j dQ.Fy{k,h,T) = FuQ , 

2,1 f ^ ^ 1 

^^ = -^ d^{k-n)F^{k,n,T) = -kF^i 

lOTT J 4 



(46) 



e^, = / an 



{k ■ h) 



1 



Fy{k,n,T) = -:^Fv2 



where we have used p^ = ZPy for the massless neutrinos. 

Integrating equations (^) and ( |40| ) over q'^dqqfo{q) and dividing them by J q'^dqqfo{q), the 
Boltzmann equation for massless neutrinos becomes 

dF^ 



dr 
dF^ 

dr 



2 ■ 4 ■ 
+ ikpFy = —-h — -{h + 6rj)P2{p) in synchronous gauge , 

+ ikpFy = 4 (0 — ikp.ip) in conformal Newtonian gauge , 



(47) 



where p = k ■ n and P2{p) = ^(3/^^ — 1) is the Legendre polynomial of degree 2. Substituting the 
Legendre expansion for Fy and using the orthonormality of the Legendre polynomials and the 
recursion relation (/ + l)Pij^i{p) = {21 + l)pPi{p) — lPi_i(p), we obtain 

Synchronous gauge — 

4 2 ■ 

Su = — Ou h , 

3 3 

Ou = k" Q^, - e,) , 



Fu2 -- 


= we, = ^e,-hF,3 + ^h + 8v 


F,i -- 


= k 





I > 3. 



(48) 
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Conformal Newtonian gauge 

Su = 

9, = 






eJ+fc^v^, 



/ 



21 



-jFu{i-i) 






l> 2. 



(49) 



This set of equations governs the evolution of the phase space distribution of massless neutrinos. 
Note that a given mode Fi is coupled only to the (/ — 1) and (/ + 1) neighboring modes. 



5.4. Massive Neutrinos 



Massive neutrinos also obey the collisionless Boltzmann equation. The evolution of the 
distribution function for massive neutrinos is, however, complicated by their nonzero mass. Prom 
equations (^), the unperturbed energy density and pressure for massive neutrinos (labeled by 
subscripts "/i" for HDM) are given by 



Ph = a 



q'^dqdnefoiq) . 



1 



q'^dqdn—fo{q) , 



(50) 



where e = e(g, r) = y/q"^ + m'^a'^, while the perturbations are 

6ph= a-^Jq^dqdnefo{q)^, 6Ph = ^a~^ J q^dqdn^fo{q)^ , 



5Tl 



a '^ Jq'^dqdnqnifo{q)^ 



'hj 



q^ 1 

q^dq dn —{uiUj - -5ij) fQ{q)^ . (51) 



Since the comoving energy-momentum relation e(g, r) depends on both the momentum and time, 
we can not simplify the calculations by integrating out the g-dependence in the distribution 



function as we did for the massless neutrinos above (see eq. 45). Instead of applying equation 
1), we expand the perturbation ^ directly in a Legendre series 



(52) 



^(fc, n, q, t) = Y^{-iy^i(k, q, T)Pi(k ■ n) . 
1=0 

Then the perturbed energy density, pressure, energy flux, and shear stress in /c-space are given by 

Ana-"^ J q'^dqefo{q)^o, 



Sph 



6Ph = -^a-^ j q^dq^Uq)-^,, 

_ - 47r _ 

{Ph + Ph)Oh = -rka~ 

{ph + Ph)Qh 



4tt _ 

T 

Svr 
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q'^dqqfoiq)^!, 
q^dq^fo{q)^2 



(53) 



17 



Following the same procedure used for the massless neutrinos, the Boltzmann equation 
becomes 



Synchronous gauge 



Wq = Wi H — h — 

3e 6 dlnq 

' e V2 -1 2 +3 + 



Conformal Newtonian gauge — 

*o = -^^i-'^-Ji > 

3e amq 

' e V2/-1 2/ + 3 +V 

Because of the g-dependence in these equations, it requires much more computing time to 
carry out the time integration for the massive neutrino. Bond k, Szalay (1983) used a 16-point 
Gauss-Legendre method to approximate the g-integration. We do not use this method and instead 
perform the integration using cubic splines (plus a remainder obtained by asymptotic expansion) 
with a g'-grid of 100 g'-points for every wavenumber k. We verified that this was enough to ensure 
a relative accuracy no worse than 10~^ by trying the integration with 200 points. Then the 
perturbations 5ph, dPh, 0^, and Qh that enter the right-hand side of the Einstein equations are 
calculated from equations (|5^) by numerically integrating ^Oi ^i and ^2 over q. 



5.5. Photons 

Photons evolve differently before and after recombination. Before recombination, photons 
and baryons are tightly coupled, interacting mainly via Thomson scattering (and the electrostatic 
coupling of electrons and ions). In Thomson scatterings, the photon energy hv is assumed to 
be much less than the electron rest mass rrie ~ 0.511 MeV and the recoil of the electron in 
the initial electron rest frame is neglected. (We are concerned with the period after neutrino 
decoupling, when T < me-) The classical differential cross section for Thomson scattering is given 
by dax/d^ = 3ot(1 + cos^0)/167r, where ax = 0.6652 x 10~^^ cm^ and is the scattering angle 
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(e.g., Jackson 1975). After recombination, the universe gradually becomes transparent to radiation 
and photons travel almost freely, although Thomson scattering continues to transfer energy and 
momentum between the photons and the matter. 

The evolution of the photon distribution function can be treated in a similar way as the 
massless neutrinos, with the exception that the collisional terms on the right-hand side of the 
Boltzmann equation are now present. We shall denote by F^(k,n,T) the momentum-averaged 



phase space density perturbation, defined as in equation (|45|) . The linearized collision operator for 
Thomson scattering is (Wilson &; Silk 1980; Dodelson & Jubas 1993) 






arieCTT 



C 



-F^{k,n,T) + F^Q{k,T) + in ■ Ve - —F^2P2ik ■ n) 



(56) 



where Ue and Ve are the proper mean density and velocity of the electrons. The term proportional 
to P2 comes from the angular dependence cos^ 9 in the Thomson cross section above. The P2 term 
was neglected by Peebles & Yu (1970), but was included by Wilson & Silk (1980, 1981). Expanding 
Fj{k,n,T) in Legendre series as in equation (|^) and using the relations h ■ Ve = i6bPi{k • n), 
Fji = kO^/A, and F^2 = 100^, the collision operator can be rewritten as 






arieCTT 



!' 



%)Pi + 9e^P2-J2(-'y^'yiPi 

l>3 



(57) 



The left-hand-side of the Boltzmann equation remains the same as for the massless neutrinos, so 
we obtain 



Synchronous gauge - 

dry 

F^2 

F-yi 



-h, 



k ( -5y — G^ ) -|- anedT 



'7/ ' 



109. 



7 



kF-y3 -\ — h + 8ri — daneO'T&'y , 



/ 



21-1 



^(/-i) 



1±1 
21 + 3 



'7(^+1) 



aneCrxF. 



7/ 5 



l> 3, 



(58) 



Conformal Newtonian gauge 



"7 

U.y 

Fy2 

F^i 



k f -d-y — Q-y] + k %!) + aneaT{ob 



^11 5 



109 

k 



21 



-kF-yS — daUeCTTQ-y , 

T^7(/-i)-^7^^7(m) 



7-3-7 

I 



aUeCTTF. 



7/ ) 



l> 3. 



(59) 



The subscripts 7 and b refer to photons and baryons respectively. 
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5.6. Baryons 

The baryons (and electrons) behave like a non-relativistic fluid described, in the absence 
of coupling to radiation, by the energy-momentum conservation equations ( ^8|) and (|29| ) with 
5P/5p = c^ = w <^ 1 and = 0. Since the baryons are very nonrelativistic after neutrino 
decoupling (the period of interest), we may neglect w and §P/5p in all terms except the acoustic 
term c'^k'^6 (which is important for sufficiently high k; note that the shear stress term k'^Q is far 
smaller so we neglect it). Before recombination, however, the coupling of the baryons and the 
photons causes a transfer of momentum (and energy) between the two components. 

From equation ( [2l| ) the momentum density T^j for a given species is related to 6 by 
ik^ST^j = (p + P)9. The momentum transfer into the photon component is represented by 
an(.aT{Ob — ^7) of equations (|58| ) and (^). Momentum conservation in Thomson scattering then 
implies that a term {Ap-y /3pb) aneCrxiO'y — 6b) has to be added to the equation for Of, (where we 
have used Pb <C Pb), so equations (p8|) and (pi) are modified to become 



Synchronous gauge 



-eb-\h, 

--e. + clk^Sb + ^anearie^-eb), (60) 

a 2,pb 



Conformal Newtonian gauge 



h = ~0b + 3<p, 

9b = --eb + cie5b + ^aneaT{e~t-eb) + k'^i^. (61) 

a 3pb 



5.7. Tight-Coupling Approximation 



Before recombination the Thomson opacity is so large that photons and baryons are tightly 
coupled, with an^o'T = t^^ ^ a/ a ^ r^^. The large size of the Thomson drag terms in the 
equations for 6^ and Of, of equations (p^)-(|6lD make these equations numerically difficult to solve. 
Therefore, in this limit we shall follow the method of Peebles & Yu (1970) to obtain an alternative 
form of the equations valid for Tc/t <C 1 and kvc <C 1. The starting point is the exact equation 



obtained by combining the second of equations (^9|) and (61) 



:i + R)eb + -Ob - clk'^5b - k^R (jd^ - ©7 ) + ^(^7 - ^fe) = (1 + ^)fcV , (62) 
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where the right-hand side is included only in the confornial Newtonian gauge; in the synchronous 
gauge it is set to zero. We have defined R = {A/3)p.y/ph. We shall see that the terms proportional 
to {9^ — 6h) and 0^ may be neglected to lowest order in max [A;Tc,Tc/r], with the result that the 
baryons and photons behave like a single coupled fluid with velocity Vb- However, we require a 
more accurate approximation to account for the slip between the photon and baryon fluids. 



From the second of equations (|59D , we have 



Ob — O-Y = Tc 



k^ -6^ - e 



7 ^7 



k'^llj 



(63) 



in the conformal Newtonian gauge; in the synchronous gauge we simply set '0 = 0. Writing 6^ as 
9b + (^7 — 9b) and using equation (|6^ , we get 

1 



db — 9-. 



l + R 



9b + k^ [ ci5b --:5y + e^] +9^- Ob 



a result that is valid in both gauges. From the third of equations (pq), we have 



e. 



Tc 

9 



+ -h + 8'n- loe^ 



-kR 



73 



(64) 



(65) 



in the synchronous gauge; in the conformal Newtonian gauge one sets h = rj = 0. We see that 



e^ ~ 



5^ X max [kvc, Tc/t] (the case kTc corresponding to acoustic oscillations with 9^ 



k5^ 



Higher moments of the photon distribution are smaller by additional powers of kTc and we shall 
neglect them in the limit of tight coupling considered here. Our goal is to obtain equations for 9b 
and 9^ that are accurate to second order in Tc- 

To get an equation for 9b we differentiate equation (U) and use equation (|6^) . Assuming that 
the gas is nearly fully ionized so that ng oc a~^ and that the baryon temperature is approximately 
the radiation temperature implying c^ ex a~^, we obtain 



Ob 



2R a 



l + Ra 



-\^b 



+ 



l + R 



-Ob --k^\=-5^ + i^]+k'\ c% - \6^ 



+ 0{tI). (66) 



This equation holds in the conformal Newtonian gauge; in the synchronous gauge one should set 
■0 = 0. Note that 6b and 6^ are to be evaluated using the flrst of equations (|5^)~(|6T|). Substituting 



equation ( pq) into equation (62) yields our desired equation of motion for max [A;rc, Tc/t] ^ 1. If 
this condition is violated, then one should use the explicit form of equations (^) and (SI) for Ob- 

To obtain an equation for 0^ we combine the explicit equations for 0^ and Ob to obtain the 



exact equation 



R-^ { 0, 



h + -Ob 
a 



cik^db + A;' 



9- 



{l + R) 



+ ^-^fcV 



(67) 



in conformal Newtonian gauge; in synchronous gauge one sets = 0. For Ob we use the 
tight-coupling approximation (substituting eq. ^ into eq. |62|) at early times and the exact explicit 



21 



equations (^) or ( plD at late times. In practice, we switch to the explicit scheme for 9b when 
Tf, = 2 X 10*^ K; we switch to the exphcit scheme for Fyi for / > 1 and Tf, = 2 x 10^ K (at earlier 
times these moments are set to zero). We have verified that these switches occur early enough to 
preserve good accuracy in the resulting photon phase space distribution. 



5.8. Recombination 

In order to compute the Thomson scattering terms in the equations of motion for photons 
and baryons we need to know the free electron density 71^(7). Our treatment is based on Peebles 
(1968; see also Peebles 1993). We summarize the procedure here. 

Because the Thomson opacity is enormous until hydrogen begins to recombine, it is sufficient 
to treat the helium as being fully neutral at all times. We define the ionization fraction of hydrogen 
as Xe = rie/riH where nn is the number density of hydrogen nuclei. The ionization rate equation 
is (Peebles 1968; Spitzer 1978) 

^ = aCr [/3(r,)(l - Xe) - nHa^^Hn)xl\ . (68) 

The factor Cr is discussed below. The collisional ionization rate from the ground state is 



Pin) = (^^^)'^%-^^/^'B^^a(2)(r,) , (69) 

where Bi = niee^ /{2ti ) = 13.6 eV is the ground state binding energy, and the recombination rate 
to excited states is 

This expression for (j)2{Tb) is a good approximation (better than one percent for T^ < 6000 K). At 
high temperatures this expression underestimates (p2 but the neutral fraction is negligible so that 
we make no significant error by setting (j)2 = for T^ > Bi/k-Q = 1.58 x 10^ K. 

Recombination directly to the ground state is inhibited by the large Lyman alpha and Lyman 
continuum opacities. Net recombination must occur either by 2-photon decay from the 2s level, 
with a rate A2s_^is = 8.227 s~^, or by the cosmological redshifting of Lyman alpha photons away 
from the line center. Peebles (1968) gives a detailed discussion of these atomic processes. The net 
recombination rate to the ground state is reduced by the fact that an atom in the n = 2 level 
may be ionized before it decays to the ground state. The reduction factor Cr is just the ratio of 
the net decay rate (including 2-photon decay and Lyman alpha production at the rate allowed by 
redshifting of photons out of the line) to the sum of the decay and ionization rates from the n = 2 
level: 

^ ^ Aq, + A2s^ls ,„.^ 

""' A„ + A2.^i, + /5(2)(r,) ' ^"^ 
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where 

/3(2)(r,)=/?(T,)e+^-'^/'=BT. ^ A„ = ^^, A„ = ^ = 1.216 X 10-5 cm , (72) 

where v^ = c/\a- For T^ <C 10^ K, it is a very good approximation to replace the number density 
Tils of hydrogen atoms in the Is state by (1 — Xe)nH- 



We integrate equation (38) using a stable and accurate semi- imp licit method with a 
large number of timesteps through recombination. Since the results are independent of the 
perturbations, we pre-compute the ionization history of a model and later use cubic splines 
interpolation to obtain Xe(r) accurately during integration of the perturbation equations. 



6. Super-Horizon-Sized Perturbations and Initial Conditions 

The evolution equations derived in the previous sections can be solved numerically once the 
initial perturbations are specified. We start the integration at early times when a given A:-mode 
is still outside the horizon, i.e., kr <^ 1 where kr is dimensionless. (We follow common usage in 
referring to waves with kr < 1 as being "outside the horizon" even though r is more appropriately 
called the comoving Hubble distance.) The behavior of the density fluctuations on scales larger 
than the horizon is gauge-dependent. The fluctuations can appear as growing modes in one 
coordinate system and as constant modes in another. As we will show in this section, this is 
exactly what occurs in the synchronous and the conformal Newtonian gauges. 

We flrst review the synchronous gauge behavior, which has already been discussed by Press & 
Vishniac (1980) and Wilson & Silk (1981), although these authors did not include neutrinos. We 
are concerned only with the radiation-dominated era since the numerical integration for all the 
/c-modes of interest will start in this era. At this early time, the massive neutrinos are relativistic, 
and the CDM and the baryons make a negligible contribution to the total energy density of the 
universe: ^totai = Pu + P-^- The expansion rate is a/ a = t"^. We can analytically extract the 
time-dependence of the metric and density perturbations h, r], 5, and 6 on super-horizon scales 



{kr <^ 1) from equations (|20|), ( [48|) and (58). The large Thomson damping terms in equations 
( pq ) drive the / > 2 moments of the photon distribution function F^i to zero. Similarly, Fi^i for 
/ > 3 can be ignored because they are smaller than Fi^2 by successive powers of kr. Equations 
(|20i| ), (|0|), (H), and (|D then give 

r^/i + Th + 6[(1 - R^,)6^ + RJ^] = , 
4 2- -1 

6, + ^e, + h = o, e,-h\6,-Ae,) = o, (73) 

15 
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where we have defined Rp = pp/{p^ + p^). For Ni, flavors of neutrinos {Ni, = 3 in the standard 
model), after electron-positron pair annihilation and before the massive neutrinos become 
nonrelativistic, Pv/p-y = (7A'',^/8)(4/ll)^'^ is a constant. 

To lowest order in fcr, the terms oc k"^ in equations ( [7^ ) can be dropped, and we have 
9i, = 9^ = 0. Then these equations can be combined into a single fourth-order equation for h: 



whose four solutions are power laws: /i oc r" with n = 0, 1, 2, and —2. From equations ( |73| ) we 
also obtain 

h = A + B{kT)-^ + C{kTf + D{kT), 

5 = {I - Rp)5^ + Rp5p = -lB{kT)-^-lc{kTf-\D{kT), 

6 6 b 

e = (i - Rp)e^ + Rjp = -lok, (75) 

o 

and A, B, C, and D are arbitrary dimensionless constants. The other metric perturbation rj can 
be found from equation ( |20a| ) : 

ij = 2C + jD{kT)-K (76) 

Press &: Vishniac (1980) derived a general expression for the time dependence of the four 
eigenmodes. They showed that of these four modes, the first two (proportional to A and B) are 
gauge modes that can be eliminated by a suitable coordinate transformation. The latter two 
modes (proportional to C and D) correspond to physical modes of density perturbations on scales 
larger than the Hubble distance in the radiation-dominated era. Both physical modes appear 
as growing modes in the synchronous gauge, but the C(A;r)^ mode dominates at later times. In 
fact, the mode proportional to D in the radiation-dominated era decays in the matter-dominated 
era (Ratra 1988). We choose our initial conditions so that only the fastest-growing physical 
mode is present (this is appropriate for perturbations created in the early universe), in which 
case 9^ = 6u = r] = to lowest order in kr. To get nonzero starting values we must use the full 
equations (|7^ ) to obtain higher order terms for these variables. To get the perturbations in the 
baryons we impose the condition of constant entropy per baryon. Using all of these inputs, we 
obtain the leading-order behavior of super-horizon-sized perturbations in the synchronous gauge: 

Synchronous gauge — 

6^ = --C{kTf, Sc = Sb = -Su = -S^, 

9^ = 0, 0^ = 0^ = _}_c(k^r^), 0, = |±i|l^^, (77) 

k^Cikrf, V-^C-^^±^^Cikrf. 
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The initial conditions for the moments ^i,l > 1, of the massive neutrino distribution can be 
related to ^o and the variables above by equations (^). To obtain the initial ^O) we can write 
the perturbed distribution function as / = fo{q){l + ^l'o) = 2/ip {exp[g/A;(T + 6T)] + 1}~^, where 
6T/T = 6u/4: by the isentropic condition. Then using equations (^^, we find the first three 
moments to be 



^0 
*2 



1^ dlnfo 
4 dlnq 



e 

qk 



din/, 







ding 



(78) 






5u 2^u 



dlnfo 
dlnq 



The higher moments ^i {I > 3) are negligible for kr <ti 1. 

The initial conditions for the isentropic perturbations in the conformal Newtonian gauge can 
be obtained either by solving equations (^), (pSJ), and (^), or using the transformations given by 
equations (|l^ and (^) (which enables us to relate the amplitudes in the two gauges). We find for 
the growing mode 



Conformal Newtonian gauge 



40C 



e„ 



V' 



15 + ARu 

9u = dc = db 

AC 
" 3(15 + 4i?^) 

20C 
15 + 4R^ ' 



-2V, 

IOC 
15 + 4R^ 



^b = -.^v 

1 



7-^7 > 



(fcr) 



-(A:^r) = -(A;V)V^, 
{krf^, 



(79) 
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The massive neutrino moments ^i in this gauge are related to by^ 9^, and 0;^ by the same 
equations (eq. ^8|) as in the synchronous gauge. As claimed earlier, iIj = (j) to zeroth order in kr 
when no neutrinos are present (i.e., R^, = 0). If we characterize the perturbations in the conformal 
Newtonian gauge by the potential ij), all matter and metric variables have a very simple form 
outside the horizon. The neutrino energy fraction R^^ enters only in the second potential as a 
result of the shear stress produced by the free-streaming neutrinos. Bardeen (1980) was concerned 
that a large shear stress would lead to large metric perturbations in the conformal Newtonian 
gauge. We see that this does not happen for isentropic growing-mode perturbations in which the 
shear stress arises solely due to the free-streaming of relativistic collisionless particles. 

We see that 5 grows with time in the synchronous gauge but remains a constant in the 
conformal Newtonian gauge before horizon crossing. Another significant difference is the larger 
value of the velocity perturbations for small kr in the conformal Newtonian gauge. Physically, this 
difference arises because velocity perturbations vanish to lowest order in the synchronous gauge 
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because the synchronous gauge spatial coordinates are Lagrangian coordinates for freely-falhng 
observers (§2). The next-order velocity perturbations differ for the neutrinos and photons because 
these two fluids have effectively different equations of state: the neutrinos are collisionless while 
the photons behave like a perfect fluid due to their strong coupling to the baryons. In the 
conformal Newtonian gauge, the lowest-order velocity perturbations do not vanish because the 
conformal Newtonian gauge spatial coordinates are Eulerian coordinates. If we were to include the 
next-order corrections to 9 proportional to A:^t^, differences between the different fluid components 
would appear in equations ([TSJ). 

In the conformal Newtonian gauge, the mode proportional to D in the synchronous gauge 
yields cj) cc if) cc 5 cc {kT)~^ . Thus, this mode corresponds to a decaying mode in the conformal 
Newtonian gauge even though it yields 8 oc {kr) in the synchronous gauge. The two gauge modes 



{A and B in eq. 75) do not exist in the conformal Newtonian gauge. 



7. Integration Results in a CDM-|-HDM Model 

We apply the results derived in the previous sections to an = 1 cosmological model consisting 
of a mixture of CDM and HDM, with parameters r^coid = 0.65, Ohot = 0.3, r^baryon = 0.05, and 
Hq = 50 km s^^ Mpc^^. The corresponding neutrino mass is rriy = 93.13 (Jlhot^^) eV = 6.985 eV. 

In Fourier space, all the k modes in the linearized Einstein, Boltzmann, and fluid equations 
evolve independently; thus the equations can be solved for one value of A; at a time. Moreover, all 
modes with the same k (the magnitude of the comoving wavevector) evolve the same way. We 
integrated the equations of motion numerically over the range 0.01 Mpc^^ < k < 100 Mpc~^ using 
41 points evenly spaced in log^^g ^ with an interval of Alog^^Q k = 0.1. The time integration was 
performed using the standard fifth- and sixth-order Runge-Kutta integrator dverk (obtained from 
netlib@research.att.com). The time integration was begun at conformal time tq = 3 x 10~^ Mpc 
with z ~ 10^ and ended at tj = 3000 Mpc with z = 13.55. The initial tq was chosen so that the 
largest k (i.e. the smallest wavelength) was well outside the horizon at the onset of the integration. 
The integration was stopped when the fluctuations were still in the linear regime, with the rms 
density fluctuation in the CDM component being about 0.2 (Ma &; Bertschinger 1994). 

The Einstein equations provide redundant equations for the evolution of the metric 
perturbations. In the synchronous gauge we chose to use ah and rj as the primary metric 
perturbation variables in the integration, and used equation ( 20b| ) and a combination of ( pOa| ) 



and (20c) as the evolution equations. In the conformal Newtonian gauge we integrated (j) using 



equation (|22b| ) and obtained ■0 algebraically using (|22d| ). In both gauges we used the time-time 



Einstein equation (eqs. |20a| and 22a) to check integration accuracy. In the conformal Newtonian 
gauge it is possible to avoid integration of the metric perturbations altogether by combining 
equations ( p2aD and ( |22b| ) into an algebraic equation for 0. However, we found that this gave 
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numerical difficulties because the initial value of (j) has to be set with exquisite precision when 
kr -^ 1. We also found it necessary to obtain the initial Os from (j) and the 5s in the combined 
constraint equations of ( p2aD and ( |22b| ) . Although the analytical expressions in equations (^) are 
good approximations for /cr ^ 1, slight deviations from the energy- momentum constraints was 
found to cause numerical difficulties. 

The photon and the massless neutrino phase space distributions were expanded in Legendre 



series (see eq. 45) with 1000 /-values in order to guarantee sufficient angular resolution. The 



massive neutrinos are computationally expensive due to the momentum-dependence in equations 



(54) and (|55|). We performed the massive neutrino calculations on a grid of 100 (/-points including 
50 /-values for every q. By setting the phase space harmonics to zero for larger values of /, we 
found spurious waves which propagate from low to high /, reflect at / = 50, and then propagate 
back to low /. These numerical artifacts are analogous to the propagation of traveling waves on 
a string with a fixed end. However, only the first three /-terms contribute directly to the source 
terms in the Einstein equations, and we are not interested in the angular power spectrum of the 
neutrinos themselves. We found that /max = 50 is adequate to ensure that the reflected waves 
from the cut-off at /max have not interfered with the low-/ harmonics. We checked that all of our 
numerical approximations are adequate by increasing the grids of /c, g, and / values as well as 
decreasing the integration timestep. We estimate that our final results have a relative accuracy 
better than 10"'^. 

We shall first present results for the metric perturbations in the conformal Newtonian 
gauge and then compare the evolution of density perturbations in our two gauges. The metric 
perturbations in the synchronous gauge have no simple physical interpretation, so we shall not 
bother presenting them. 

Figure 1 shows the time evolution of the metric perturbations (f){k,T) and 'ip{k,T) in 
the conformal Newtonian gauge for all 41 values of k. The overall normalization was chosen 



arbitrarily (corresponding to C = —1/6 in eqs. 77 and U^). The difference between tp and (j) in 



the radiation-dominated era is due to the shear stress contributed by the relativistic neutrinos 
(eq. 79) which make up a fraction Ri, = 0.4052 of the total energy density. On scales much 



smaller than the horizon, ip corresponds to the Newtonian gravitational potential and (p = ip 
in the matter-dominated era. As is well known, the potential is constant for growing-mode 
density perturbations of CDM in an Einstein-de Sitter universe. In a mixed dark matter model, 
however, the CDM density perturbation growth can be suppressed by the lack of growth of HDM 
perturbations so that tp decays slowly. 

The behavior of the metric perturbations can be understood as follows. All of the 41 
A;-modes are outside the horizon at early times when r < 0.01 Mpc. The horizon eventually 
"catches up" and a given /c-mode crosses inside the horizon when kr is about tt. The modes 
with larger k (i.e., shorter wavelengths) enter the horizon earlier. If a given /c-mode enters the 
horizon during the radiation-dominated era, the tight coupling between photons and baryons due 
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to Thomson scattering induces damped acoustic oscillations in the conformal Newtonian gauge 
metric perturbations, which are exhibited in Figure 1 by the modes with k > 0.1 Mpc^^. (In 
fact, it is not the speed-of-light horizon that sets the scale for the oscillation and damping of the 
potential. Rather, as we show below, it is the acoustic horizon. These horizons are similar during 
the radiation-dominated era because the sound speed of the photon-baryon fluid is c/Vs.) The 
modes with k < 0.1 Mpc~^ enter the horizon during the matter-dominated era and do not oscillate 
acoustically because the Jeans wavenumber fcj = {AirGpa'^ /cf)^''^ has then become much larger 
than the wavenumbers under investigation. 

We can understand the oscillations more quantitatively by studying the Einstein equations 
(I24 ) in the conformal Newtonian gauge. Analytical solutions can be found for a perfect fluid with 
no shear stress, in which case (j) = tp. Using c^ = 6P/6p = p/p and a/ a = 2t~^ /{I + 3c^) (from 
eqs. 18 and ^), equations (22) can be combined to yield 



^'^ + ttS^^'^ ^ ^^"'^^''^ " ° ' ^^°^ 



whose solutions are Bessel functions with a power-law pre-factor: 



/>± = (fcc,T)-'^J±,(fec,r), ^^ 2(l + 3c2) - (81) 



The (/>_ solution corresponds to the decaying mode discussed in §6, so we shall ignore it. In the 

I and z^ = f , 



radiation-dominated era, c^ = o and z^ = |, so that 



,0/9 , , f constant , /cCcT ^ 1 , , , 

4., = ikc.r)-^''j,M>^c.r) « { ^_, ^^^^^^_^^ ^^_^ ^^ ^ ; (82) 

This type of behavior is apparent in Figure 1 for r < Teq. This analytic solution holds of course 
only in the absence of neutrinos; obtaining the correct amplitudes for ip and (p in CDM-I-HDM 
models shown in Figure 1 required the full integration discussed in this paper. 

After the universe becomes matter-dominated, acoustic damping of the potential ceases, 
and the only physical process causing the potential to change is the free-streaming damping of 
perturbations in the massive neutrinos. We shall discuss this further after examining the evolution 
of the density perturbations. 

Figure 2 shows the evolution of the density perturbations for the five particle species in the 
two gauges from our numerical integration. Three wavenumbers are plotted: k = 0.01 Mpc~^ 
(Fig. 2a), k = 0.1 Mpc~^ (Fig- 2b), and k = 1.0 Mpc~^ (Fig- 2c). Each mode is normalized with 
the same initial amplitude for (p as in Figure 1. There are several notable features: 

Before horizon crossing — 

(1) The initial amplitudes of the (5's are related by the isentropic initial conditions: 
6^ = 6u = 5h = 4:6h/3 = 45c/3- The behavior of 5 outside the horizon is strongly gauge-dependent. 
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In the synchronous gauge, Figure 2 shows that ah the (5's before horizon crossing in the 
radiation-dominated era grow as a? . This confirms equations ( [77| ) since aij^ oc r at this time. It 
is straightforward to show that in the matter-dominated era (for 0=1), 5 oc r^ oc a for all modes 
before horizon crossing. In the conformal Newtonian gauge, the (5's remain constant outside the 
horizon as derived in equations ([79[). After horizon crossing, however, the perturbations come 
into causal contact and become nearly independent of the coordinate choices. As one can see, 
(5c, (5f,, and bh in the two gauges are almost identical at late times. In fact, one can show using 
equations (^) and (pqa) that if -(/; = constant, then (5(Syn) — 5(Con) = 2^, accounting for the 
slight differences apparent in the figures. 

After horizon crossing — 

(2) For CDM, the /c-modes that enter the horizon during the radiation-dominated era 
behave very differently from those entering in the matter-dominated era. The critical scale 
separating the two is the horizon distance at the epoch of radiation-matter equality (oeq ~ 10~^): 
/ceq = 27r/req ~ 0.1 Mpc~^ for our parameters. For the modes with k > 0.1, horizon crossing 
occurs when the energy density of the universe is dominated by radiation; thus the fluctuations 
in the CDM can not grow appreciably during this time. For the photons and the baryons, the 
important scale is the horizon size at recombination (orec ~ 10^'^): fcroc = 27r/Trec ~ 0.025 Mpc~"^. 
The modes with k > 0.025 (see Figs. 2b and 2c) enter the horizon before recombination, so the 
photons (long-dashed curves) and baryons (dash-dotted curves) oscillate acoustically while they 
are coupled by Thomson scattering. The coupling is not perfect. The friction of the photons 
dragging against the baryons leads to Silk damping (Silk 1968), which is prominent in Figure 2c 
at a ~ 10~^'^. The baryons decouple from the photons at recombination and then fall very quickly 
into the potential wells formed around the CDM, resulting in the rapid growth of 5b in Figures 2b 
and 2c. 

(3) Neutrinos decouple from other species at T ~ 1 MeV and a ~ 10"^''. At this early time, 
both the massless and the 7 eV massive neutrinos behave like relativistic coUisionless particles. 
The massive neutrinos become non-relativistic when SksTi, ~ rui,, corresponding to Om ~ 7 x 10~^ 
for 7 eV neutrinos. Close inspection of the figures at a ^ Onr reveals that 6fi (short-dashed curve) is 
indeed making a gradual transition from the upper line for the radiation fields to the lower line for 
the matter fields. Although the Jeans length of a fluid is not well defined for coUisionless particles 
such as the neutrinos, the criterion for free-streaming damping is similar to the Jeans criterion 
for gravitational stability: free-streaming is important for k > kf^, where A;|,(a) = AnGpa^ /{v'^) 
and {v^)^'"^ is the neutrino velocity dispersion. When the neutrinos are relativistic, t> ~ 1 and 
kis{a) oc a~^ for a < Ocq. After the neutrinos become non-relativistic, the neutrino speed is given 
by (i;2)i/2 ^ s/cBT^/m^ = 3kBTo^u/am,y , implying (t;2)i/2 ^ l5a~i(m^/10eV)-^ km s~^ In the 
matter-dominated era, we have AnGpa^ = ^HqU'^ from the Friedmann equation, and therefore 

kUa) = 8a'/^ l^^^hMpc-K (83) 

In Figure 2a, since horizon crossing occurs when the free-streaming effect is already unimportant. 
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the evolution of 5h is very similar to that of CDM. In Figures 2b and 2c, however, the free 
streaming effect is evident and the growth of 5^ is suppressed until A;fs(a) grows to ~ k. After 
kis{a) > k, 5h can grow again and catch up to 6c- Since /cfg oc a^'^, the larger k modes suffer more 
free-streaming damping and 5h can not grow until later times. The damping in 6h also affects the 
growth of Sc, slowing it down more for larger f^hot compared to the pure CDM model. This effect 
is also apparent in the power spectra shown in, for example. Ma & Bertschinger (1994). 



8. Summary 

Physical quantities are independent of the coordinate systems they are computed in. For 
historical reasons, most calculations of linear fluctuation growth have been carried out in the 
synchronous gauge. In this paper we explored an alternative gauge, the conformal Newtonian 
gauge, which is free of the gauge ambiguities and coordinate singularities associated with the 
synchronous gauge. We derived the coordinate transformation relating the two gauges and 
presented the linear theory of isentropic scalar gravitational perturbations in parallel for both 
gauges. The complete set of evolution equations are given: the Einstein equations for the metric 
perturbations, the Boltzmann equations for the photon and neutrino phase space distributions, 
and the fluid equations for CDM and baryons. 

The use of the conformal Newtonian gauge was motivated by our work on the HDM initial 
conditions in CDM+HDM models (Ma &: Bertschinger 1994). In order to sample the neutrino 
phase space accurately for an N-body simulation, we followed individual neutrino trajectories 
by numerically integrating the geodesic equations in a perturbed Friedmann-Robertson- Walker 
background metric. The conformal Newtonian gauge proved to be the most convenient choice 
for this calculation because the geodesic equations have simple forms and the coordinates do not 
become seriously deformed at late times. 

In this paper we applied the linear theory to the CDM+HDM model under study: r^cow = 0.65, 
f^hot = 0.3, Jlbaryon = 0.05, and Hq=50 km s~^ Mpc^^. The evolution of the density fields for all 
five particle species in both gauges was presented. We also illustrated the gauge dependence of the 
density fields before horizon crossing and discussed the physical interpretation of the results. 

Interested users may obtain our programs to integrate the perturbation equations by 
contacting one of the authors. 
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Fig. 1. — The scalar metric perturbations (f){k,T) (Fig. la) and ip{k,T) (Fig. lb) in the conformal 
Newtonian gauge as a function of r. The 41 curves from left to right correspond to 41 values 
of k between 100.0 Mpc~^ and 0.01 Mpc~^. The labels Tnr, Teq and Tree indicate the time 7 eV 
neutrinos become non-relativistic, the matter-radiation equality time, and the recombination time, 
respectively. 
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Fig. 2. — Evolution of the density fields in the synchronous gauge (top panels) and the conformal 
Newtonian gauge (bottom panels) for 3 wavenumbers k= 0.01 (Fig. 2a), 0.1 (Fig. 2b) and 1.0 
(Fig. 2c) Mpc~^. In each figure, the five lines represent 6c, Sb, S^, by and bh for the CDM (solid curve), 
baryon (dash-dotted), photon (long-dashed), massless neutrino (dotted), and massive neutrino 
(short-dashed) components, respectively. 



